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We use harmonically trapped systems to find the leading gradient corrections of the superfluid local 
density approximation (slda) - a density functional theory (dft) describing the unitary Fermi gas 
(ufg). We find the leading order correction to be negative, and predict the coefficient of the 
long-range static response = 1 .5(3) - a factor of two smaller than predicted by mean-field theory - 
thereby establishing a new and experimentally measurably universal constant. 
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UNIVERSALLY describing two-component Fermi sys- 
tems with short-range interactions of infinite scat- 
- tering length Ug — )• 00, the unitary Fermi gas 
(ufg) [1] not only approximates the dilute neutron mat- 
ter found in neutron stars [2], but is directly realized in 
cold-atom systems [3], allowing experiments to bench- 
mark many-body techniques used to study astrophysical 
phenomenology. Despite the simplicity of the system - 
the lack of length-scales, for example, implies that the 
equation of state E (n+ ) oc n^^^ - the system is strongly 
interacting and admits no known perturbative expan- 
sions: A quantitative description requires experiments 
or ab initio computations. 

Ab initio techniques, however, can only address a few 
questions - direct quantum Monte Carlo (qmc) simu- 
lations, for example, can study systems with at most 
a few hundred particles. It is therefore imperative to 
benchmark computationally tractable models of macro- 
scopic phenomena so that they can be used to answer 
outstanding phenomenological questions, such as the 
origin of glitching in neutron stars [4] . 

Density ftinctional theory (dft) is an in principle exact 
approach, widely used in nuclear physics (see [5] for a 
review), and in quantum chemistry to describe normal 
(i.e., non-superfluid) systems. It provides a framework 
capable of assimilating ab initio and experimental re- 
sults into a computationally tractable and predictive 
framework. In this letter, we extend one such dft - the 
superfluid local density approximation (slda) - to de- 
scribe the inhomogeneous behaviour of harmonically 
trapped systems. We use dft to analyse recent experi- 
mental and theoretical results, noting discrepancies and 
the asymptotic behaviour toward the thermodjmamic 
limit, and establish the leading order gradient correc- 
tions to the slda which we find to be negative. The slda 
then uniquely predicts the low-energy static response of 
the UFG to quadratic order, which is crucial for a proper 
low-energy description of the ufg: this therefore makes 
significant progress towards a predictive framework for 
studying superfluid phenomenology. 

At low-energies, the ufg can be characterized by a 
superfluid effective field theory (eft) describing phonon 



dynamics [6J. The eft admits a controlled power- 
counting scheme: The leading order (lo) contains a 
single dimensionless parameter - the Bertsch param- 
eter [7] £, which characterizes the equation of state 
£(n) = £,£fg ('"•+) where £fg = S/Sn+Ep is the energy 
density of a free Fermi gas with the same total density 
n+ = kp/(37t-^), and Ep = h^kp/2Ta is the Fermi energy. 
At next-to-leading order (nlo), two additional dimen- 
sionless coefficients {6, ^8J appear which characterize the 
static and d5mamic low-frequency and low-momentum 
response. We shall address only the value of the static re- 
sponse Xq as defined by adding small external potential 
6Vq(x) = 5cos(qx) to a homogeneous system: 



6n+(x) =Xq5Vq(x) + 0{62). 
To NLO in the superfluid eft, the response is 
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where is a universal dimensionless constant fg]. This 
normalization for is numerically close to unity (c^ = 1 
for non-interacting fermions), appears simply in the 
energy of trapped fermions (see Eq. ||6}), and is inde- 
pendent of the £, and pairing parameters in the slda. 
The other universal constant Cuj describes low-energy 
dynamical properties, and enters through the phonon 
dispersion relation fio] a)q/(qcs) = 1 + Ca)q^/(24f,kp) -|- 
0(q^ In q) where Cj = hkj: is the speed of sound. 

While many techniques have been employed to cal- 
culate the Bertsch parameter £, w 0.37 (see [11] for a 
survey), nlo coefficients have only been considered in 
a few cases: The e-expansion [12] (expanding in spa- 
tial dimension: e = 4 — d) gives « 8/5 -I- O(e^) and 
Cx ~ Co) + O(e^) [ [131, while Bogoliubov-de Gennes {sdc) 
mean-field theory [8) finds = 7/3 and Coi = 0.7539. 

The EFT breaks down for small systems and near the 
boundary of clouds, so to connect with finite-size qmc 
calculations, we use dft. The Hohenberg-Kohn theo- 
rem [14I asserts the existence of a universal functional of 
the density alone whose minimum describes the ground 
state of the ufg. The exact form of this functional is 
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non-local and unknown, but a local formulation - an 
extended Thomas-Fermi (etf) functional [Ql [i5|fi8ll - 
describes well some energetic and dynamical aspects of 
the UFG. It fails, however, to properly describe finite-size 
effects in homogeneous systems [ [TgllSol . 

To describe these properties we use the Kohn-Sham 
formulation [|2ij which includes an auxiliary kinetic 
density t+. While this is formally equivalent to the 
Hohenberg-Kohn formulation, the addition of a kinetic 
density allows a local formulation to describe finite-size 
features of the system. In particular, the finite-size prop- 
erties of non-interacting systems are exactly reproduced. 
Interacting versions have been considered [ST, 23], but 
one finds that the finite-size effects are not properly sup- 
pressed IjigJ l2oil . The suppression can be realized by 
include an additional auxiliary anomalous density, v, 
representing the pairing field Jig) I20I [24] , resulting in 
the slda: 
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Here a is the inverse effective mass, |3 is the self -energy, 
y controls the pairing, and 5A characterizes the leading 
order gradient term (known as a Weizsacker correction). 
The unitary limit is realized when we take the wave- 
vector cutoff kc — > 00 to infinity (see Ref. P4] for details). 
In homogeneous systems, the gradient corrections van- 
ish, and one can use the equations in the thermodynamic 
limit to replace the parameters |3 and y by the more phys- 
ically relevant quantities f,s and r\ = A/Ep, where A is 
the pairing gap (see the appendix of [20] for details). 
When applied to inhomogeneous systems, however, one 
must hold the parameters (3 and y fixed to define the 
functional. 

This form ||3| subsumes earlier dfts. In particular, the 
well-studied Bdc mean-field equations are reproduced 
with unit effective mass a = no Hartree 

term |3 = 0, and no gradient corrections 5A = 0. The 
Kohn-Sham form discussed in |22I |231 neglects the v = 0, 
while the etf form |T2l [T5I is reproduced if one neglects 
both the anomalous density y and the kinetic density 
T+. As discussed in [19], none of these restricted forms 
can even qualitatively characterize the finite-size effects, 
but we still consider the etf functional as it is much eas- 
ier to solve numerically while retaining the asymptotic 
properties of trapped systems: 



h^[]/4+ 5A) (Vn4 
8m n+ 
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The leading gradient term here derives from a semi- 
classical expansion of the kinetic energy [25] with an 
additional Weizsacker correction 6A. Superfluid hy- 
drodynamic phenomenology [17] and vortex dynam- 
ics (appendix |E| suggest that 6A = 0. The resulting 



ETF is completely determined by the value of £,. A 
simple calculation [.18J shows that the etf model has 

= Co, = 9/4 + 9 5A = 9/4. 

The SLDA was originally constrained by qmc calcula- 
tions of the continuum state, and validated with qmc 
calculations in a harmonic trap I26I [27} . These valida- 
tions, however, provided only a weak test of the slda 
form. In particular, the symmetric thermodynamic limit 
does not provide enough information to constrain the 
effective mass, and the original variational trap results 
were not sufficiently accurate to exhibit the appropriate 
scaling in the thermod5mamic limit [ [TgllSol . 

Recently, experimental and ab initio qmc results for 
homogeneous matter in the continuum and in periodic 
boxes were used to more rigorously test the form of 
the SLDA [20"]: the best fit to current unbiased results 
is consistent with £, = 0.3742(5), a = 1.104(8), and ti = 
0.651 (9). Here we estimate the leading order gradient 
correction 6A by reconsidering the energies of trapped 
systems. 

The static response in the thermodynamic limit can be 
calculated using the same techniques as in the Bdc [8J 
and one finds 1. 28 1 



= ja-|-9 5A, 



(5) 



which is independent of £, and r|. This demonstrates 
how the effective mass and gradient corrections play a 
similar role, as pointed out in p9)[3o[ . 

From the nlo superfluid eft |'6'^'8l, one finds the 
energy of the ufg in an isotropic harmonic trap with 
trapping frequency w to depend on the coefficients £, 
and 
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hcu 
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This form naturally suggests the abscissa x = (3N + )^^/^ 
so that the asymptotic behavior of E is linear: we prefer 
to use the square of the energy E^, 



16E^ 



h2a,2(3N + )8/3 
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as f, appears on the intercept, and appears directly. It 
is interesting that, in the non-interacting system, shell- 
effects appear at the same linear order x, leading to a 
fundamental uncertainty in the coefficient f ^ ~ 
1.7. Pairing suppresses these shell effects, yielding a 
well-defined asymptotic slope c^, which can also be 
determined from the semi-classical approximation (see 
appendix [B|. 

In Fig. IiTwe display qmc results for trapped unitary 
systems. The dotted lines guide the eye through sev- 
eral variational bounds obtained using fixed-node qmc 
(fnqmc) calculations. In these methods, one avoids any 
sign problem by sampling a restricted set of wavefunc- 
tions with the same nodal structure as an initial reference 
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Figure i. (color online) Various qmc and experimental results 
for trapped systems (with lines to guide the eye). The results 
with dotted lines are from fnqmc calculations that provide 
only upper bounds due to the nodal constraint. From top to 
bottom these results are from |26 | (blue), [27] (red), and [20] 
(green). The (magenta) points with solid lines are from a 
lattice calculation [31] that is in principle unbiased. The large 
(black) crosses for N+ e {4,6} are from [32]. The solid light 
(grey) curve shows the shell structure for free fermions in 
the trap (the curve has been shifted down from 1 to facilitate 
comparison). Finally, we include the latest experimental value 
for £, = 0.376(5) from [33] as a (yellow) band and the best fit 
value of £, = 0.3742(5) to all homogeneous ab initio qmc results 
from [20] at X = 0. Coordinates have been scaled as in (|6} to 
demonstrate the scaling. The corresponding particle numbers 
N+ are listed along the top axis, and emphasize the closed 
shells which occur for N+ e {2,8,20,40,70}. (Note that all 
methods agree for the point N+ = 2 (not shown) which admits 
an analytic solution.) 

ansatz. By improving the ansatz and varying the pa- 
rameters, these bounds have come down over time, and 
the lowest (green) curve Hoi represents the best boimd 
to date. Note that this is the only set of results that 
demonstrates the expected linear scaling ||6| predicted 
by the effective theory. We suspect that numerical issues 
or the nature of the ansatz in the other cases introduced 
spurious lengths scales that violate this scaling (see [20] 
for further discussion.) 

The solid (magenta) line guides the eye through cal- 
culations based on lattice techniques [31]. In principle, 
these are unbiased ab initio results, but it is somewhat 
troubling that most lie significantly above the variational 
bounds. They also display large shell effects that are 
virtually absent in the fnqmc results. For comparison, 
we have included the energies of free particles shown 
in Fig. [4] as a light (grey) curve, shifted down from 1 
to facilitate comparison. As we shall see, although the 
non-interacting slda reproduces these shell-effects, the 
interacting slda exhibits a marked lack of shell effects, 
consistent with the fnqmc. The lattice results thus seem 
qualitatively inconsistent with the others. 

A third variational technique [ [32] based on a corre- 
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Figure 2. (color online) Dft models compared with data 
from Fig. [1] The thick solid (black) curve is the Slda with 
parameters f, = 0.3742(5), a = 1.104(8), and n = 0.651 (9) but no 
gradient corrections 6A = 0. The thin dashed (blue) curve is 
the ETF with the same £, = 0.3724 and 6A = 0. For comparison, 
the upper thin solid (red) curve corresponds to the best fit etf 
model (£, = 0.468, 6A = —0.164) described in fi"5l to the fnqmc 
result [27] . The linear asymptotic forms ij = E, + c-^x with 
Cx = 7a/3 2.58 (slda) and = 2.25 (etf) are shown as short 
thin lines extending from x = 0. The inset is a log-log plot of 
the deviation the etf model makes from the asymptotic form 

y = £, + 9/4x + ax^/^ H where a f« —2.2: the solid curves are 

the deviations (y — above (blue) and —{y — £, — 9/4x) below 
(black)) while the dotted lines are the expected order (9/4x 
above and qx^/^ below). 

lated Gaussian approach provides very tight bounds, but 
is limited to small systems. These are shown as (black) 
crosses for N + € {2, 4, 6). Unfortunately, at these three 
points, all methods agree, and significant discrepancies 
between the lattice and fnqmc results only appear at 
larger N + . 

Finally, we include the latest experimental results as a 
(yellow) band [33] and the best fit value of f, = 0.3742(5) 
to all homogeneous ab initio qmc results from [20] at 
X = 0. According to the scaling ||6}, the results should 
approach this point. 

In figure [2] we overlay the dft results for the trapped 
systems. The best fit slda without gradient corrections 
(solid) and the etf with the same £, = 0.3742 and 5A = 
(dashed) have almost exactly the same structure, demon- 
strating the ability of the kinetic term in the slda to 
model the gradient effects of the etf. While the dft 
results appear to approach a linear asymptotic form, the 
actual slope determined from the static response is 
only realized for extremely large N + . This is expected 
due to the very slight suppression of higher order terms 
~ x^/^: to suppress these corrections by an order of mag- 
nitude requires 0.1 w x^/^ which implies that N+ ^ 10^. 
Thus, there is virtually no hope of directly extracting 
the value of from the qmc simulations presented in 
Fig. [1] without using a model to extrapolate to large N + . 
Both dft models ultimately exhibit this behaviour as we 
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Figure 3. (color online) Slda model as in Fig. [2] with various 
values of the gradient corrections A added perturbatively. 



demonstrate with the etf in the inset. 

Note that the slda predicts higher energies for sys- 
tems with small particle numbers - especially the N + 6 
{2,4, 6} systems where all ab-initio methods agree. To cor- 
rect for this, the leading gradient correction needs to be 
negative A < 0. We explore these effects in Fig. Rlby per- 
turbatively including the gradient for various values of A. 
(Note: these corrections are less than 15% for all systems 
as shown in Fig. [5] The remaining corrections from a 
fully self -consistent solution will not significantly alter 
these results.) One could also increase the effective mass, 
but to match even the modest correction of A ^ —0.1 
requires TUgff ^ 1.4m. which spoils the description of 
homogeneous systems [[ig} goj and quasi-p article dis- 
persions [34]. Higher order gradient corrections might 
help, however, there are several different gradient cor- 
rections - each requiring additional coefficients - and 
insufficient ab-initio results to constrain these. These 
neglected terms will not affect the coefficient c^- Fi- 
nally, there is the possibility that the functional could 
be generalized as discussed in ["Tg"], but the success of 
the three-parameter slda [20] suggests that corrections 
along this line would be small. 

We note that a negative gradient correction is some- 
what surprising since a naive expansion of an attractive 
non-local interaction — V(x — y)n(x)n(y) yields a positive 
gradient correction (see appendix |F}. 

Thus, we conclude from Fig. [3] that the slda will re- 
quire A w —0.12(3) to describe both homogeneous boxes 
and trapped systems. The slda therefore predicts 



1.5(3) 



(7) 



where the error is approximate (not to be taken as a stan- 
dard deviation). This is about half the value predicted 
by the Bdc mean-field calculation, though the effect in 
the static response oc c^/E,, is cancelled by the ex- 
cessively large sdc value of f, = 0.5906. Intuitively, c^/f, 
is four times larger than in the non-interacting system. 



indicating that momentimi-dependent density fluctua- 
tions are suppressed by interaction, though this effect is 
mostly due to the reduced value of £,. 

As discussed, cannot be directly extracted from the 
QMC results without a model capable of extrapolating 
well into the thermodynamic limit. To extract more 
directly, one should consider systems that minimize the 
sensitivity to the breakdown of the eft at the boundary 
of the system. To do this, consider how the density n = 
n(x) depends on a smoothly varying potential V = V(x): 



^ _ Cx (VV)2+4(H-V)V^Vh.^ 

m 



64 (^-V)3 

,3/2 



n-TFlxj 



3^2 (2£,) 3/2 



m\3/2 



(8) 



This is valid to nlo in a static background with constant 
phase (i.e. not near a vortex). Thus, one can directly 
search for deviations from the Thomas-Fermi (tf) profile 
that are sensitive to gradients in the potential. Applying 
a modulation V oc cos(qx) directly probes the static re- 
sponse and is feasible in qmc simulations. Optical latices 
could be used similarly in experiments, however, mea- 
suring the density to sufficient accuracy is likely to be a 
challenge - to compensate for the numerical suppression 
0^/64 while avoiding contamination from higher-order 
terms will require percent level accuracy. Experiments 
should thus probably retain traps with axial or spherical 
symmetry so they can benefit from averaging techniques 
like the inverse Abel transform to reduce noise in n(x). 
Adding a small dimple to the core of the trap with vary- 
ing widths will allow experiments to probe the response 
while retaining the benefits of averaging to reduce noise. 

Confirming the value of will provide an important 
validation of the slda fimctional, and provides another 
benchmark for models of the ufg: To reliably predict 
low -energy behaviour of the ufg, a model should repro- 
duce the LO and nlo coefficients - £,, c^, and Cw - of the 
superfluid eft 
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Appendix A: Contradictory QMC Results 

As noted in Fig. [ij there are unresolved contradictions in 
the QMC results for larger N+ with the lattice results fyi] 
exceeding the fnqmc variational bounds [20]. It would 
be useful to have an alternative method calculate the 
energies for N+ = 8 and N+ = 12 to resolve between 
these: 

^ = 11.64;0]06 ^I"^ 12.01(2) (N+=8) 
ncu u.iz^ 

^ = 20.765;°-°45 %3Mc ^ (N+ = 12) 

Appendix B: Semiclassical Expasion 

The superfluid eft is closely related to the well-studied 
semiclassical expansion (in h.) [25] for non-interacting 
systems, and one can derive similar expressions to those 
arising from the eft. The utility of the eft is to organize 
the universal coefficients for the interacting superfluid 
system. 

N+ 
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Figure 4. (color online) Shell effects in the trapped non- 
interacting two-component gas. The scaling is the same as 
in Fig. [1] The asymptotic bounds | ^ C;^ £ 1 .7 have been 
drawn as dotted lines. 

The shell effects for non-interacting systems is shown 
in Fig. [4] (this was shifted down to match £, for compari- 
son in Fig.jiJ. The semiclassical expansion systematically 
organizes contributions from volume effects, surface ef- 
fects, periodic orbits, etc. It thus provides some insight 
into the breakdown of the eft: shell for the harmonic 
trap, for example, effects appear at the same order as the 
NLO corrections. These are shown in Fig. I4] demonstrat- 
ing a fundamental uncertainty in the slope 5 ^ c-^ ^ 1 .7. 
Considering the static response allows one to extract the 
non-rnteractrng value of = 1, which lies in this band. 
The fact that these corrections appear at the same order 
is related to nearby breakdown of the eft formula, and 
the subsequent slow approach to the asymptotic scaling 
in the thermodynamic limit. What is non-trivial is that 



pairing acts to sufficiently suppress these shell effects 
so that a well defined slope emerges in harmonically 
trapped systems. Perhaps this can be explained within 
the semiclassical theory, but the author is not aware of 
such a discussion in the literature. 

Appendix C: Perturbative Gradient Corrections 

In Fig. [3] we show the size of the gradient corrections 
Egrad as a percent of the total energy E for the 6 A = —0.15 
required to bring the slda in rough accordance with 
the smallest trapped systems. This demonstrates that 
the gradient corrections may be included perturbatively, 
simplifying the numerical calculations. 
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Figure 5. With Egrad/^ (ir* percent) for 6A = —0.15, demon- 
strating that the correction is indeed perturbative. 

In order to perform a fully self-consistent calculation 
of the gradient corrections, one must modify the single- 
particle self-energy to include a term 

^'""^ ~ 8m \n+ ) ~ &m y ^ n+ 

Accurately computing the derivatives - especially in the 
tails of the cloud - presents a mild numerical challenge, 
and will only alter the energies at the percent level, and 
so is not required for the present analysis. This self- 
consistency will be important, however, if one wants 
to compare the density profile of the smallest trapped 
systems with qmc results. 

Appendix D: Asymptotic Behaviour 

Figure |6] is a full-sized version of the inset of Fig. [2] This 
shows on a log-log scale the deviation of the numerically 
computed y (x) from the etf model from the expected 
asymptotic form ||6| y = f, -I- c^x — 2.2x^/^ -I- ■ ■ ■ . (We do 
not claim anything universal about the coefficient 2.2 
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Figure 6. (color online) Asymptotic behaviour of the etf model 

with £, = 0.3742 and 6A = 0: ij = £, + 9x/4 + qx^/s h where 

a !=i —2.2. We plot the deviations y — f, (upper solid blue curve) 
and — (y — £, — 9x/4) (lower solid black curve) along with the 
next terms 9x/4 (upper dotted blue line) and l.liJ^^ (lower 
dotted black line) to demonstrate the asymptotic powers. 



here. It was simply obtained numerically by fitting the 
asymptotic form: likely, a more careful analysis of the 
SLDA results will yield a different value.) This confirms 
that the etf does approach the asymptotic form, but 
only for N+ ~ 10^ as expected from the proximity of 
the correction. We also show the same analysis for the 
SLDA in Fig. [7 however, since it is much more expen- 
sive numerically, we have not approached the N+ w 10^ 
threshold (note the different scales). Although the agree- 
ment does not seem as good, this is simply a function of 
the small particle numbers: a comparison with the etf 
results in the same region shows similar deviations. 



Appendix E: No Weisacker Term in the ETF Model 

Superfluid hydrodynamic simulations with the etf find 
phenomenologically that the etf without gradient cor- 
rections works best to describe collisional dynamics in 
the UFG [17]. Here we also argue that they must van- 
ish to give a sensible description of vortices. The etf 
model Q is equivalent to a modified Gross-Pitaevskii 
Equation (gpe) [16 35] with a complex field ¥ describ- 
ing dimers of mass 2m normalized such that the density 



i9t¥: 



4m 



+ 2(^EF(n+)-^i]- 6A 



m|¥| 



Consider a single stationary vortex V oc e^'t'f (r) embed- 
ded in a uniform gas with background chemical poten- 
tial |j. = £,Ep(n.;^) = bf^/^(r = 00). The phase yields a 
centrifugal term: 



^l±l^^l!^ + bfV3(,) 

4m 



H f(r)=0. (El) 
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Figure 7. (color online) Asymptotic behaviour of the slda 
model with a = 1.104(8), £, = 0.3742(5), and t| = 0.651(1): y = 

£,+ 7ax/3 + ax^/^ H where a ^ —1.1. We plot the deviations 

y — £, (upper solid blue curve) and — (y — £, — 7ax/3) (lower 
solid black curve) along with the missing correction terms 
7«x/3 (upper dotted blue line) and 1.1x7^^ (lower dotted black 
line) to demonstrate the asymptotic powers. Note that the 
scale is quite different from figure |6] owing to the additional 
computational complexity of the slda. The approach to the 
asymptotic form is consistent. 



The Weizsacker term modifies the effective mass in the 
gradient term, but does not similarly alter the centrifugal 
term (the first term in Eq. | |Ei[ |) since it acts only on the 
modulus 1^1 . The vortex develops a cusp near the core: 
let f (r) = ar" where awl; then 



Qr«-2 1 + 4 5A 2 a-2 



4m 



4m 



^^7/3^7a/3 _ ^^^0 



The terms on the left-hand-side are divergent and must 
cancel, fixing = 1 -I- 2 6A. With no Weizsacker term 
6A = 0, this yields the familiar a = 1, but the presence 
of a Weizsacker spoils the cancellation between the gra- 
dient and centrifugal terms, and the density profile of 
a vortex at the core becomes non-analytic with a cusp 
of a fractional power. This non-analytic cusp causes 
unphysical dynamical evolution of the vortex. 



Appendix F: Naive Gradient Corrections 

To show that a negative gradient correction is somewhat 
surprising, consider expanding a local interaction in 
terms of the separation f = (x — y )/2 in the spirit of the 
density-matrix expansion (dme) 



V(x 



-y)n(x)n(y) 
V(2f)n^(R)- 



- V(2f)n(R + f)n(R-f) 
2V(2f)[f- Vn(R)]^ 



Thus, a naive expansion for an attractive potential would 
imply a positive gradient correction. It is apparent that 
this simplistic argument does not apply to the strongly 
interacting ufg. 
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